Vortex States and Phase Diagram of Multi-component Ginzburg-Landau Theory with Competing 

Repulsive and Attractive Vortex Interactions 
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We investigate the behavior of vortices of multi-component superconductivity, realized in MgB 2 and Fe-based 
superconductors, within the framework of Ginzburg-Landau (GL) theory in terms of numerical calculations of 
the time-dependent GL equations and the variational method. It is revealed that close to the critical point of the 
composite system the inter-component coupling makes the system behave as a single component superconduc- 
tivity in most cases. However, when the bare mean-field critical points of the two components coincide with 
each other, and furthermore the inter-band coupling disappears at the same temperature, interesting phenomena 
occur as follows. Vortices interact attractively at large separation and repulsively at short distance in certain 
parameter space. Because of the non-monotonic interaction profile, phase separations between vortex clusters 
of triangular order and the Meissner state take place, which indicates a first-order phase transition associated 
with the penetration of the magnetic field into a superconductor sample. Phase diagrams of vortex states are then 
constructed with the associated magnetization curve. It is found that all these behavior interpolates the features 
of the type I and II superconductors. 

PACS numbers: 74.25.Ha, 11.27.+d, 74.25.Uv, 74.70.Ad 



I. INTRODUCTION 

A quantized vortex is a topological excitation of superflu- 
idity and superconductivity as the hallmark of a single wave 
function for the macroscopic quantum state. Revealed first 
by Abrikosov based on the celebrated Ginzburg-Landau (GL) 
theory, a vortex in superconductor carries a normal core with 
size of the order of superconducting coherence length £, and 
a quantized magnetic flux distributed over the London pen- 
etration depth A; in superconductors categorized as type II 
with k - A/i; > 1/ V2, as opposed to type I with k < 1/ v2, 
vortices repel each other because of the circular supercurrents 
and thus form a triangular vortex lattice. In type I supercon- 
ductors, vortices would attract each other in order to gain the 
superconducting condensation energy, and thus collapse into 
a continuum of normal region at equilibrium. 

The recently discovered multi-band superconductors, such 
as MgB2 1 1 ] and ironpnictide superconductors! 2 1, can exhibit 
novel phenomenaOHD ( see Ref. lfTOl for a review), without 
counterpart in single-band superconductors. It was first pro- 
posed theoretically by Babaev and Speight[ 11, 12 1 that, when 
the London penetration depth falls in between the two coher- 
ence lengths, vortices may attract at large separation and re- 
pel at short distance because of the competition of different 
length scales in different condensates. It was then discussed 
by Moshchalkov et al. that this situation was realized in their 
sample of MgB2, as manifested by the unconventional stripe 
and gossamer-like vortex patterns| 13 1. They coined the term 
of type 1.5 superconductor for the class of multi-band super- 
conductors. 

While this scenario is intriguing and has attracted consider- 
able interests ifHUTrjl . there are several points waiting for fur- 
ther investigation. In the theoretical side, one notices that the 
original proposal was based on weak inter-component cou- 
pling limit, and the discussions were implicitly extended to 
low-temperature regime. In doing so, one should always keep 
in mind that GL theory works only close to the critical point 



of the composite system, at least in principle. Any result ap- 
pearing only away from the critical point needs a careful and 
independent check. In the experimental side, the observed 
inter-vortex separation ~ 2//m is much larger than the distance 
associated with the possible energy minimum estimated by the 
theory [llj, lacking a consistent understanding. The random 
vortex configurations reported in Refs. l"L3l[T7l seem incom- 
patible with common experiences in material science: parti- 
cles governed by the Lennard- Jones interactions form solids 
with regular orders, and thus those random patterns are hard 
to be considered as a property of equilibrium at usual condi- 
tions in absence of random pins. 

In the present work, we reveal first a case that the bare 
mean-field critical points of the two components coincide with 
each other, and furthermore the inter-band coupling disap- 
pears at the same temperature. We introduce two vortices 
into a square superconductor by setting appropriate boundary 
condition, and find with the approach of time-dependent GL 
(TDGL) that, for appropriate parameters in the GL free energy 
functional, the two vortices are separated by a distance asso- 
ciated with an energy minimum, which does not change with 
the system size for simulations. This verifies the attraction 
between vortices at large distance and repulsion at short dis- 
tance. The vortices are found stable against thermal fluctua- 
tions. We evaluate the full dependence of interaction potential 
on vortex separation by the variational method. Distribution 
of vortices are then simulated based on the Langevin dynam- 
ics. Instead of uniform vortex lattice in type II superconduc- 
tors as well as the lamella structure in type I superconduc- 
tors, we observe phase separations among clusters of triangu- 
larly ordered vortices and Meissner regions, which indicates a 
first-order phase transition associated with penetration of vor- 
tices into the system upon increase of external magnetic field. 
Based on these observations we construct the phase diagram 
for superconductors with the novel vortex interaction. 

We also investigate the vortex interaction when two compo- 
nents have different bare critical points and a finite inter-band 



coupling. Close to the critical point of the composite sys- 
tem, the superconductivity in different components is strongly 
correlated and only one divergent length associated with the 
variation of order parameters at the vortex core exists. Thus 
vortices are purely repulsive or attractive close to the critical 
point. 

It is noticed that novel vortex attractions were discussed 
40 years ago in single-band superconductors, which caused 
phase separation between Meissner phase and vortex lattice, 
and discontinuous jumps in magnetization, especially at low 
temperature (see Ref.[18| for review). These superconduc- 
tors are characterized by small GL numbers, and the attractive 
vortex interactions were attributed to correction to the GL ap- 
proach from the BCS theory] 19 20|, which is different from 
the system discussed by flTT] [T2) and addressed in the present 
work. 

The remaining part of the paper is organized as follows. In 
Sec. II, we present the free energy functional and the structure 
of a vortex. In Sec. Ill, we present the results obtained by 
numerical simulations of TDGL. In Sec. IV, the interaction 
potential is obtained by the variational method. In Sec. V, 
the vortex configuration is simulated based on the Langevin 
dynamics. Sec. VI is devoted discussions, and the paper is 
concluded by Sec. VII. 



We introduce the dimensionless quantities for convenience, 

x = A lX ', % = <¥ w % A = AiH lc V2A', T = %T' , 

y = Y\ ai \, B = 7/ lc V2B', J=^|fJ', 

(2) 

where *Pj = \a\ \ jfi x is the bulk value, H lc - J An \a\\ *Pf the 
thermodynamic critical field of the first condensate. B is the 
magnetic induction and J is the supercurrent. 

For clarity, we drop the prime in the dimensionless quan- 
tities in the following calculations. Then we have the free 
energy in the dimensionless units 



r = E i=w gjroa+ft^ + al^v-AjW 



+(VxA) 2 - 7 (W 2 +W l ), 



(3) 



where k\ - Ai/gi is the GL parameter. 

Minimizing T with respect to ^F, and A, we obtain the GL 
equations 



T, \ IT, | -'4', i ( — V- A 



¥i - y x ¥ 2 = 0, 



(4) 



- S * 2 + £ l*2l 2 ¥ 2 + ^(^V-AW 2 -yF 1 = o, (5) 

a x Pi m 2 \iK l ) 



Vx VxA 



II. MODEL 



The GL free energy functional with interband scattering is 
given by 



r = z,- = 



1.2 



a,ra 2 +fi^i 4 + ^|( 



-mv 



2e 



A H 



+i(VxA) 2 -r(^ 2 + l P*'P 1 ), 

(i) 

with all symbols conventionally defined ll2T1 . Here we fo- 
cus on the Josephson-like inter-band coupling with strength 
y, noticing that the main results remain qualitatively the 
same even including other possible interactions. For MgB2 
y > 0, while for iron-based superconductors, y is presumably 
negative El . The temperature dependence of the quadratic 
terms are aff) = a,(0)(l - T/T c ) with a,(0) < 0, and 
y(T) = y(0)(l - T/T c ). We notice that in this case the physics 
is the same for the whole temperature regime starting from T c 
except a renormalization of length. We leave the discussion 
on more general cases to Sec. VI. 

For convenience, we define the following lengths /I, = 
-\JmjC 2 fii/87r\aj\e 2 and ^,- = y/H 2 /2mj\aj\, the penetration depth 
and coherence lengths in the respective single-band conden- 
sates (y = 0). For the present interest, we adopt the co- 
efficients in the GL free energy functional Eq. (fT) such 
that at T - 0, %\ - 51nm, A\ - 25nm, £2 = 8nm, and 
Ai - 30nm at very low temperatures. The inter-band coupling 
isy(0) = -0.4ai(0). 



2lK[ 



(vj/*vT' 1 -TjWD-PFirA 



I {^2^2 - ^W*) - |Y 2 | 2 A) . 



(6) 



Equation (J6l) describes the screening of the magnetic field by 
the superconducting condensates. Using the London approx- 
imation, the effective London penetration depth for two-band 
superconductors is 



•I" I/Wl^l0l 2 + — 1^201 

ni2 



(7) 



where ¥,0 is the bulk value of the /th superconducting con- 
densate. The response of two-band superconductors to mag- 
netic fields is described by a single length scale A, because 
both condensates couple to the same gauge field. For y - 0, 
we have Tjo = 1 and ¥20 = ya^i/ci/^- The interband cou- 
pling shifts the bulk value and thus modifies the corresponding 
penetration depth. 

Next we construct a vortex solution to Eqs. (HI), p) and 
Q, It is observed that the two condensates should have the 
same vorticity and phase around the vortex core to save en- 
ergy if y > 0, while for y < the phase shift should be n. For 
condensates with different winding number, the vortex is frac- 
tional quantized and the energy diverges logarithmically with 
vortex size in bulk|4|, and thus thermodynamically unstable. 
Without loss of generality, we focus on the positive y > as 
in the case of MgB2[ 1 1. Presuming a straight vortex line, and 
we look for a vortex with the following structure 



% = fi{r)e™ and A = 



na(r) 
Kir 



-e e , 



(8) 



where r is the distance from the vortex core, eg the unit vector 
in the azimuthal direction and n the vorticity. Substituting the 
ansatz into Eqs. (HI), pi and d6]l, we have 

- Mr) + ffr) - - 2 fe/i + -a r fi) + n {a ~ l) h - yfi = 0, 

(9) 



a ±Mr) + %fHr) 



^^f 2 + id r f 2 ) + '^f 2 )- 7 f^0, 



dta - -d,a + \f} + — fi I (1 - a) = 0. 
r \ m 2 ' 



(10) 



(11) 



In the limit of r — > oo, the wave functions recover the bulk 
values f w and f 20 . Defining f 20 = 77/10 with 77 > 0, we have 
the equations for /o and 77 



-l+/, 2 o-W = 0, 



a 2 #2 3,-, . „ 

— 7+^-»7 (i +y?)-7 = o. 



(12) 



(13) 



The radial variation of the wave functions and vector poten- 
tial in the asymptotic region of r — > 00 can be found and are 
given by 



/ = a/1 +7V + c f 1 ex P - 



V%J 



IjSi /ar 2 r, 
■/2 = VF — + _ + co exp 



(14) 



(15) 
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FIG. 1: (color online). Left: schematic view of two vortices in a su- 
perconducting square disk. Right: numerical discretization scheme, 
f'" is denned on nodes, J, A and U are defined on bonds, and B- is 
defined inside the plaquette. 



III. TDGL APPROACH 

To study the interaction between vortices, the structure of 
the nonlinear vortex core is important and numerical calcula- 
tions are necessary. We first calculate the structure of vortex 
and their interaction by minimizing the free energy in terms 
of TDGL equations defined in the following way IflOl |23l 



h 2 2e 6T 
(d, + i—tyVi = — — + C u 

InijDi ' ft ' ' 6*¥* t( 



-(-d,A + VO) 
c c 



5T_ 
'<JA 



+u, 



(19) 



(20) 



with D, the diffusion constant, <x the normal conductivity, and 
<t> the electric potential. By choosing a proper gauge, we can 
take <£ = 0. C,j represents thermal noises satisfying \lj\ — 

and (£/(xi,?i)f/x 2 , h)) = I>5(xi - x 2 )6(h - t 2 \ with 7= 1,2 
and A. Since we are primarily interested in the equilibrium 
properties, detailed dynamic relaxation process is irrelevant. 
Here, time is in units of t\ - £;\ID\, the normal conductivity 
<x in units of c 2 Ti/4nA^. 



1 + c a exp 



A v 



(16) 



At large distance, there is only one length scale for the two 
condensates. The penetration depth can be obtained straight- 
forwardly 



/ \m 2 p 2 \«i 77/ 



(17) 



To calculate the coherence length, we substitute Eqs. ( 14 1 and 



( 15 1 into Eqs. (J9l) and ( 10 1 and linearize the equations. £,, is 
equivalent to the length scale of the small fluctuations in the 
bulk, and is given by the largest solution to the equation 



2 + 3777 - 



1 



2*?#J 



ay y m\ 1 
2— +3-- ' 



77 m 2 2 K \g) 



-y 2 =0. (18) 



It is clear that the interband scattering modifies the coherence 
length and penetration depth in a nontrivial way. 



A. Numerical techniques 

In order to integrate the TDGL equations, we general- 
ize the numerical method developed for single component 
superconductors [24 1 to two-component superconductors and 
solve the TDGL equations Eqs. ( fT9] ) and (|20j> numerically. For 
the parameters we are interested, the line tension of vortices is 
high, therefore we approximate vortices as straight lines. The 
problem is then simplified into two dimensions. To maintain 
the gauge invariance after discretization l25l . we use the link 
variable defined as 



UJx,y) = exp 






wci AJ&dtu 



fin 



(21) 



where fi is either x or y. Then the TDGL equations can be 
rewritten as 

cW~i Yj ^^(^*i)-^i+l*il 2 *i-r*2-^i =0, (22) 

1 n=x,y 



quette with the index (i, j) is 



- I p=x,y 



<xd,A + Vx VxA = J + &, 



(23) 



(24) 



+t (d^WtW - t/^ 2 ^(f/;^)]) . - ' 

The disk is partitioned into square meshes of size h as 
schematically shown in Fig. [T] The magnetic field in the pla- 



B,,i,j 



l — W ■ ■ 

1^ With W z , y = f/; ; iJ+1 C/;. yf/. v; w c/, ; i+hj , 

(26) 



and the supercurrent 






(2W2> 



i ™i /ff mWu/W _ 77* »t/U)itM>* \ 



(27) 



where t 1 * denotes the £th condensate. The y component is 
obtained similarly. After the discretization, the TDGL equa- 
tions become 



yd) 



c(l)|2WD 



(2) 



r(l) 



dX> = (i - iwpytf 1 + ypv 1 + C + 



t-U*. 



w- 



(28) 



m%Dl 



I^^S = (" 



Iml 



A np(2),2Nvp(2) , ,Ap(l) -I- / (2) -I- "" ' 



\j +U ', ,-ijl 



-2>F W 



jT.-' +£/* . . ,*-' , 



-2¥ u 



/l 2 



A 2 



(29) 



' *. *»7 a - ■*' ( >7 ^ /i 2 A > *.7 jj (+1J m 2 ■*' '■■•J i,j i+ij ^i.j J ' 



(30) 



d t U y;iJ = -±U y . u 1m[- 



»'-. 



>• 'J ij i.i+i 






(^)l 



(31) 



Equations ( 28 1, ( 29 1, ( 30 1 and ( 3 1 1 are integrated by the for- 
ward Euler method. 

Vortices are introduced through the boundary condition 



A(r + Li) = A(r) + V X i, T*(r + L,) = <F*(r) exp(/27^,/O ), 

(32) 
with / = x, y and ^ x — H a Ly and ^- v = 0. H a is the applied 
magnetic field and should obey the vortex quantization con- 
dition via §d\ ■ A = 2m<t>o with an integer m, which yields 
H a = 2m<t>o/L 2 . Here Oq = /ic/2<? is the flux quantum. 



B. Vortex structure 



We minimize T in Eq. (1) by solving the TDGL equa- 



tions numerically under the boundary conditions ( 32 1 speci 



fying the number of vortices in the system. The structure of 
a single vortex obtained by the TDGL equations for m — 1 
is shown in Fig. [2ja). Three characteristic length scales are 
evident. 



C. Vortex attraction 

We then introduce two vortices into a square disk with size 
10/li < L < 4QA\. The disk is large enough for one to ne- 
glect the finite-size effect. We find a solution with two vor- 
tices separated by r m - 2.7 A\, independent on L, as shown in 
Fig.[2jb). This indicates unambiguously an attractive interac- 
tion at large separation and a repulsive interaction at short dis- 
tance between vortices. The minimal energy associated with 
this vortex separation is therefore demonstrated to be an in- 
trinsic property of the superconductor. We also confirm the 
stability of the above vortex solution against thermal fluctua- 
tions in the present TDGL approach. 

The size of the magnetic flux lies between the radii of the 
normal cores associated with the two order parameters. As 
two vortices approaching each other, they attract first by over- 
lapping their normal cores associated with *Fi as shown in 
Fig. [2ja). When the two vortices get closer, strong repulsion 
caused by the magnetic flux sets in. An equilibrium config- 
uration is reached by compromising the attraction and repul- 
sion, where the normal core associated with *P; is overlapping 
strongly while ^2 is not as shown in Fig.[2jb). Thus vortices 
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FIG. 3: (color online). Vortex structure obtained by the TDGL equa- 
tions (symbols) and the variational calculations (lines). 



attract at large separation and repel at short distance. It is the 
competition between the sizes of the normal cores associated 
with the two components and the area of magnetic flux accom- 
panying the vortex cores that governs the interaction between 
vortices. 



IV. VARIATIONAL CALCULATIONS 

To obtain the full dependence of the interaction energy on 
the vortex separation, one needs to fix vortices at desirable 
positions. Since vortex is of structure and not a point object, 
evaluation of vortex interaction at a given distance beyond the 
London limit is not straightforward [ 27 1 . In order to overcome 
this difficulty we resort to the variational method developed by 
Jacobs and Rebbi |28 1 and generalize it to two-band supercon- 
ductors. The essence of this approach is to construct "good" 
trial functions for two vortices at given separation. 

A. Vortex structure 

To reproduce the asymptotic behavior of a vortex at r — > oo, 
we use the following trial functions for one vortex 



/t(r)= Vl^W + exp -- — 2( /l -' r 7 /! )' (33) 



FIG. 2: (color online), (a) Profiles of order parameters, magnetic 
field and supercurrent for a single vortex, (b) A stable two-vortex 
solution in a finite sample derived numerically using the approach of 
TDGL equations. The vortex separation is r m = 2.1 A\. 



\M«i 7/ I V%Jl=o 



)J( /2 ,,r'//!),(34) 



(r)=l+exp(--MJ( fl , ///!), (35) 



where f\j, fcj. and ai are variational parameters. In the limit of 
r — * 0, fi = and a — > r 2 according to Eqs. d9b, ( flOp and (Hi. 

We have f lfi = - VTT^, / 2;0 = - ^f (i + 9' fl() = -1 
and «i = aolA v . For a giant vortex with vorticity equal to 2, 
we have / u = / 1-0 /( V2£,) and / 2 ,i = / 2 ,o/( V2£ v ) in addition. 
Other coefficients are variational parameters to be determined 
numerically. 

Denote the whole set of variation parameters f\ j, / 2 /, a/ and 
so on by V. The GL free energy is of fourth order in V,- and 
has the form of 



-id 



(2), 



r(3) T 



r = T + £,■ T^'Vi + Efej nf V ' V J + Sfejs* *?* *W* 

(36) 
We use the Newton method of optimization ll29l with the iter- 
ation procedure 



y(m+\) _ y(m) _ V~< fjj 11 . jj(m) 



(37) 



where the superscript (m) represents the value at the mth 
step. Dj - dT /dy ; \ v=v m, and the Hessian matrix H t j = 

d 2 T I dV idV {\ v =v t»i). The stationary solution to Eq. (B7li cor- 
responds to the (local) minimum of the free energy. 

Following the procedure of Eq. ( 37 1, we obtain the vari- 



ational coefficients, from which we can construct the vortex 
solution. We truncate the higher-order corrections to the trial 
functions at n — 6 and find the solution of a single vortex with 
vorticity 1. The results reproduce well those obtained by the 
direct minimization of the GL free energy functional as shown 
in Fig. [3] 



B. Vortex interaction 

We proceed to consider the interaction between two vor- 
tices. For convenience, we introduce complex representation 
of the 2D coordinates z — x + iy, and then the winding phase 
part in Eq. (J8J) becomes exp(/0) = -\Jz/z* ■ To construct trial 
functions for two-vortex solution, we follow Ref. [28 1 and use 
the conformal transformation of the complex plane 



z = (z'f - (d/2) 2 



(38) 



It is straightforward to see that the origin in z plane have two 
images in the z' plane at ±d/2, and that when z' varies by 
2jt, z varies by An. This means that we may map a one-vortex 
solution in z' to a two-vortex solution in the z plane. The phase 
factor of the two-vortex solution with vortices cores at z — 
±d/2 can be constructed by this transform. We seek the wave 
functions of the form 



Vi(z,z') 



1/2 



fi(z,z*), (39) 



with i = 1,2. The trial function for fi can be constructed by 
the following consideration. For d — > °°, two vortices behave 



independently, while at d ~ two vortices merge and form 
one giant vortex with vorticity 2. In addition, we need also 
a term to describe the interaction between two vortices. The 
trial function therefore can be constructed 



Mz,z*) = a>f 2 (\z-i\)jf»(\z + j\) 

+{l-u>ym^tf ) (\z\) + 6f i (z,?), 



(40) 



where f ' and /. ' are single-vortex solution with vorticity 1 
and 2 respectively obtained by the variational calculations in 
the previous section, and 6fi accounts for the interaction. u> 
interpolates two independent vortices solution and one giant 
vortex solution. The factor in the second term at the right- 



hand-side of Eq. ( 40 1 is to ensure that the wave function van- 
ishes at the vortex cores z - ±d/2. The interaction contribu- 
tion may be constructed as follows 



Sfi(z,zl 



cosh(V2^|z|) 



ztozu/w^fey+fiy 



(41) 

The first factor is again to make sure that wave function van- 
ishes at the vortex cores, and the second factor accounts for 
the fact that the interaction vanishes when z — > °°. 

Using the similar reasoning parallel to the construction of 
the trial functions of *!*,, we can obtain the trial function for A 



A = w[ 



Ki(?-d/2) 



(1) 



+3a-oV 2 >(izi). 



(Ml) 



Kite»+d/2) 



(1) (k + f 



' (z, z*) , 



(42) 



where a m and a (2) are for the single-vortex solutions with vor- 
ticity 1 and 2. The interaction contribution has the form 



6a(z,z*) 



1 



cosh(|z|) 



[zai (z,z*) + z* a 2 (z, z*)] ■ 



with 



a k (z,z*) - ^^«£, 

i=0 7=0 



|Z| 

iJ ~2 



2/ 



:?H7 



(43) 



(44) 



with k - 1,2. 

Once the trial functions are constructed, we can do the op- 
timization according to Eq. (37i. We put two vortices at 



{-d/2, 0) and {d/2, 0), and then calculate the distribution of 
the wave functions, magnetic field and supercurrent. As dis- 
played in Fig. HI when two vortices approach from d - 8/li, 
the condensate in the first band overlaps first, which causes 
attraction between vortices, see Figs. |4](a) and (b). As d is 
reduced further, the magnetic fields start to overlap and strong 
repulsion sets in as shown in Fig. |4| c )- Finally, two vortices 
coalesce into a giant vortex with vorticity 2, see Fig. |4](d). 

The resultant separation dependence of interaction between 
two vortices is displayed in Fig. B] The agreement between 
the estimates on the position of minimal energy derived by the 
variational technique and by the TDGL equations (Fig. |2|b)) 
serves a successful check of the validity of the variational cal- 
culations. At large separation r, the attraction decreases ex- 
ponentially and the saturated energy corresponds to the self 
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FIG. 4: (color online). Distribution of the order parameters (purple and blue lines) and magnetic field (yellow line) at vortex separation (a) 
d = 8/li, (b) d = 6-ii, (c) d = 3/li, and (d) d - 0/ti. Distribution of the supercurrent at vortex separation (e) d = %X\, (f) d = 3/li, and (g) d = 0. 



energies of two isolated vortices. For r < r m , the repulsion in- 
creases sharply and the energy at r - corresponds to a giant 
vortex with vorticity 2. 
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FIG. 5: (color online). Dependence of the total energy E per unit 
length on the separation r between two vortices. The distance r„, cor- 
responds to the energy minimum of the two-vortex system. Another 
length denoted by r is the lattice constant when the ordered vortex 
patterns are formed as discussed later. 



C. Comparison on numerical approaches 

To calculate interaction between vortices poses a chal- 
lenge to theory since vortices are extended objects. In the 
London limit, the interaction of vortices has been calculated 
analytically [21]. For general cases, one has to introduce con- 
straints to fix two vortices at a desired separation. A legitimate 
procedure is to fix vortices through the boundary condition 
Eq. ( |32"| ). Vortices in a disk of type II superconductor tend to 
keep away from one another, but they cannot leave the disk 
as imposed by the boundary condition. By minimizing the 
free energy, we obtain the interaction energy and the equilib- 
rium vortex separation for a given system size L. Gradually 
increasing L, we obtain the interaction energy versus the vor- 
tices separation. However, for vortices with attraction at large 
distance and repulsion at short distance, this approach cannot 
give the dependence of the interaction energy on vortex sepa- 
ration since the distance between two vortices is always fixed, 
corresponding to the minimum of the interaction energy. We 
notice that this approach is still the most legitimate way to 
prove the existence of attraction at large distance and repul- 
sion at short distance. 

In numerical calculations of the GL free energy, one might 
alternatively introduce pinning to vortices by fixing the am- 
plitude and/or phase of superconductivity order parameter in 
a certain region near the vortex cores [30]. However, this 
method may introduce artifacts when two vortices are close to 
each other since the the order parameters change dramatically 
near the vortex core. Furthermore, the local constraints are 
sometime insufficient to pin vortices when the interaction be- 
comes strong when vortices are close. In contrast, in the varia- 
tional approach shown above the global structure (such as core 
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of the vortices and asymptotic behavior far from the core) of 
the vortices is maintained, and one only adjusts the detailed 
structure through the variational calculations. If the trial func- 
tions are appropriately chosen, the variational approach gives 
superior results. 



V. PHASE TRANSITION AND PHASE DIAGRAM OF 
VORTEX STATES 

A. Vortex configuration 

The novel interaction profile shown in Fig. E\ is expected to 
dominate the vortex state of a macroscopic system, and thus 
the phase diagram. In order to elucidate the situation, we per- 
form numerical simulations using the overdamped Langevin 
dynamics |31| 



rjdr/dt = -dE/dr + F 



(») 



(45) 



where E is the pair potential in Fig. pi and F (n) is the white 
noise force with (F") = and {Ff\i)F ( "\t')) = 2Tt}8 uj 6{t - 
f) with T being the temperature and 77 the viscosity. 

Initially vortices are randomly distributed corresponding to 
a high temperature and the system is annealed to T - by 
gradually decreasing T, which yields the ground state. Equa- 



tion ( 45 1 is solved by the 2nd order Runge-Kutta method and 



the simulation box is divided into cells with a cutoff radius r c 
to accelerate the simulation. We fix the number of vortices N v 
and set the simulation box with aspect ratio 2 : V3 to accom- 
modate the triangular lattice. The magnetic induction is given 
by B - 2N v Q>o/ V3^ 2 with L being the length of the simula- 
tion box. Periodic boundary condition are used and dt - 0.02, 
r c - 7.8 and N v - 400. The results are checked successfully 
by using finer dt, larger cutoff and more vortices. 

A typical ground-state vortex configuration for small aver- 
age magnetic induction is presented in Fig. RJa). A circular 
droplet of vortices with internal triangular order is obtained. 
This vortex configuration is clearly determined by the vortex 
interaction. In the presence of attraction, vortices prefer to 
stay together to form a cluster. For a vortex at the cluster sur- 
face, the number of neighbors is smaller than that inside the 
cluster, and thus it bares a higher energy. This results in a 
positive surface tension for the vortex cluster, similar to type 
I superconductors. The circular cluster minimizes the surface 
energy. On the other hand, the repulsion force at short dis- 
tance prevents the vortices from collapsing. Due to the repul- 
sive force, vortices inside the cluster are triangularly ordered, 
same as type II superconductors. The circular cluster of the 
closest packing triangular lattice minimizes the total free en- 
ergy. The lattice constant tq is slightly smaller than r m (see 
Fig.Bl) due to the contributions from vortices at large separa- 
tions where interaction is attractive. 

Besides the circular droplet of vortices, vortex stripe and 
vortex void are also observed at intermediate densities which 
minimize the free energy for the system with given (finite) 
size and number of vortices. At a small vortex density, the 
droplet phase [Fig. |6|a)] with triangular order is stable. Upon 



increasing the magnetic field, the vortex droplet expands, and 
at a threshold field a vortex stripe phase [Fig. |6|b)] is pre- 
ferred. The stripe then evolves into vortex void configuration 
[Fig. |6lc)] when B is increased further. When B becomes 
even larger, the void structure disappears and a perfect trian- 
gular lattice [Fig. |6|d)] emerges and remain stable until the 
superconductivity is broken completely at H C 2- 

Since the stable vortex configuration presumes the minimal 
surface area, the transition fields between two configurations 
can be evaluated by comparing the surface areas associated 
with the configurations. For a vortex droplet, the surface en- 
ergy is 2nRcr s with the radius of the cluster R — J *\/3N/2k 
and the surface tension cr s . For a stripe configuration, the sur- 
face energy is 2Lcr s . The energy consideration gives the tran- 
sition field at B^ - (bo/nr^. Similar argument gives the fields 
of other structure transitions: transition from vortex stripe to 
vortex void at B sv = (2n V3 - 3)<I>o/(3;rcTy), and transition 

from vortex void to vortex lattice at B v \ - 2<J?o/ V3r2. The nu- 
merical results are consistent with these analytical estimates. 
It is noted that these transition fields depend on the shape of 
the simulation box. In the thermodynamic limit, the circular 
vortices droplet is the only ground state at low magnetic field. 



B. Phase transition and phase diagram 

The phase separation between the Meissner region and 
a vortex cluster indicates unambiguously a first-order phase 
transition, at which clusters of vortices penetrate into the sys- 
tem upon increasing the external magnetic field [11|. The 
transition magnetic field H c \ is given by H c \ - <\n(e - 
n\ A)/Oq, with n\ - 3, e and A defined in Fig.|5]for the energy 
of a single vortex line and the energy drop associated with the 
vortex attraction per unit length. The contribution from the 
attraction component is small since e » A (see FigBl, and 
therefore H c \ is close to that of a type II superconductor with 
flux-line energy per unit length e. For large magnetic fields a 
uniform vortex lattice of triangle (Fig(6fd)) has been observed 
same as type II superconductors. 

Now we can construct a mean-field phase diagram of vor- 
tex states in the two-component superconductors. The H - T 
phase diagram is depicted in Figj7ja), where the upper critical 
field H C 2 is the same of a type II superconductor determined 
solely by the condensate with the smaller coherence length. 
The B — T phase diagram is given in FiglTlb), with the bound- 
ary between the phase separation and the uniform triangular 
vortex lattice given by B c \ - «2 c £ > o/'"m w i m n 2 — 2/ V3. 

It is illuminating to compare the response to applied mag- 
netic field in the present system with those in conventional 
type I or type II superconductors. For type II superconduc- 
tors, an extremely dilute vortex lattice penetrates into a su- 
perconducting sample at H c \ and the magnetization curve is 
continuous at the penetration. For type I superconductors, the 
magnetic field penetrates into the sample and breaks the su- 
perconductivity at H c associated with a jump in the magneti- 
zation curve. For superconductors with competing repulsive 
and attractive inter-vortex interaction, clusters of vortex lat- 
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(a)B = 0-046Q> Q fA{ 



(b)B=0.094O /^ 2 
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FIG. 6: (color online). Vortex configurations at several typical values of magnetic field: (a) vortex droplet, (b) vortex stripe, (c) vortex void, 
and (d) triangular vortex lattice. Black dots denote the centers of vortices and open region is the superconducting region. Red lines are Voronoi 
polygons. The size of simulation box is (a) L = lOCUi, (b) L = 70/li, (c) L = 50/1], and (d) L = 42/1!, with the number of vortices N v = 400 
fixed. 



tice with given lattice constant ro penetrate into the sample 
associated with a discontinuous jump in magnetic induction 
from zero to B c \. It then increases gradually with the external 
magnetic field until H c i at which the superconductivity is de- 
stroyed (see inset of FigJ7|». Therefore, the magnetic behavior 
of these superconductors interpolates those of the type I and 
type II superconductors. 



C. Effect of thermal fluctuations 

It has been revealed that thermal fluctuations should be 
weak in MgB2|23]. In the present case, the competition be- 
tween long-range attraction and short-range repulsion may en- 
hance thermal fluctuations, and warrants additional treatment. 
Phase diagrams of particles with a hard-core repulsion plus an 
attractive tail have been investigated intensively. It is known 
that the first-order gas-lattice phase boundary starts from T = 
and p - 0, and that all other transitions, such as gas-liquid, 
liquid-lattice and possible lattice-lattice transformations, take 
place at temperature k^T jE a ~ O(l), where E a is the strength 
of attractive potential (see for example |32|). In the present 




FIG. 7: (color online). Mean-field H -T (a) and B - T (b) phase di- 
agrams of superconductors with competing inter-vortex interaction. 
The red (thin) lines indicate the first-order phase transition, and the 
blue (thick) lines are for the second-order phase transition. Inset is 
for the dependence of magnetic induction on the applied field. 



flux-line system, the typical flux segment relevant to distor- 
tion of vortex lattice is estimated as L z =* ~\Jea 2 /A with a =* r m , 
which minimizes the energy associated with the tilt modulus 
and vortex interaction (see for example [33]). This gives an 
effective attraction potential e e g - L Z A - VeAa 2 . To esti- 
mate e and A, we neglect the inter-band coupling (y - 0) for 
simplicity. The attraction is caused by the overlap of the con- 
densate *P] , and thus has an order of the energy of normal core 

^ ~ (It) i? - ^he energy per unit length of a single vortex is 
contributed from the normal cores of two condensates and the 
magnetic energy, and is given bye^(^) ji + (^) m (;§)■ 
Assembling these results and taking into account the mean- 
field temperature dependence for the lengths, the condition 
ksT - e e g for various phase transitions to occur [32 1 spells 
as (1 - T/T c ) =* fGj, where / ~ 100 for a system with 
coherence length and penetration depth of the same order, 

Gi = ?( t^ fl ri»)i I is m e Ginzburg index for the first compo- 
nent with %\ (0) the coherence length at zero temperature and 

H lc = JAnc^lBi. Since G, ~ 10~ 6 as shown in Ref. ED, the 
thermal fluctuations are weak except the very small regime 
close to T c , i.e. 1 - T/T c ~ 10~ 4 , where gas-liquid, liquid- 
lattice and lattice-lattice transformations may be possible. 



Effects of thermal fluctuations have also been investigated 
by Langevin dynamics. In order to avoid possible artifacts 
due to insufficient annealing, we heat the system from the 
ground state with vortex configurations such as that shown 
in Fig.[6ta). At finite temperatures the perimeter of the cluster 
wiggles while both the cluster itself and the inner triangular 
order remain stable. Only at temperatures close T c , single vor- 
tices are evaporated from the cluster surface. These simulation 
results are consistent with the above estimate on thermal fluc- 
tuations in the present system. Through all the simulations, 
we cannot find any random vortex patterns, such as gossamer- 
or stripe -like ones lfT3l . 
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VI. DISCUSSIONS 

Here we discuss a more general case that each condensate 
has individual mean-field critical point T cl in absence of cou- 
pling and a finite inter-band coupling. The inter-band scatter- 
ing couples two bands and enhances the critical temperature 
to T c > T ci . JH 

For T < T ci and weak inter-band coupling y «c 1, the super- 
conductivity is realized independently, and the physics is qual- 
itatively the same as that of two decoupled bands y - O. irTTl 
For T c \ < T < T C 2, the superconductivity in the band 1 is in- 
duced by the the band 2 through the Josephson coupling. It 
is the regime discussed in Ref. Ifl2l . When T C j < T < T c , 
the superconductivity is purely induced by the inter-band cou- 
pling. The physics can be quite different in different tempera- 
ture regimes. We notice that in principle the GL theory is valid 
only close to T c , and thus in the highest temperature region. 

From Eq. (TTJ, T c is given by the condition a\{T)a2(J) - 
y 2 - 0, where a,(T) depends on the temperature and can be 
derived from the BCS theory (35j. Close to T c , r = ^3p = 
(T - T c )/T c <c 1, the order parameters up to the leading order 
can be derived from Eqs. (HI and pi for a uniform supercon- 
ductor [141 



< = 



a 2 a\T 



a\Pi 



afc' 



¥ 



d\azT 



20 



a\fi\ + alfc ' 



(46) 



(47) 



We then consider some deviations from the bulk values W, = 
¥,0 + (pi. From Eqs. Q and Q, the deviations are governed 

by 



QT101 + 3p{¥ 



10<Pl 



2m i 



7<p2 = 0, 



(48) 
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FIG. 8: (color online). Distribution of the order parameters and mag- 
netic field at several typical temperatures. The results are obtained 
with T cl = 0.6, T c2 = 0.8 and T c = 0.87. 



a 2 (f>2 + 3^2^20^ - ^—V 2 4>2 - 74>\ 
2m 2 



(49) 
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Taking the solution of form <p-, ~ exp[-x/( V2£ v )], we can de- 
rive the coherence length which is divergent at r — > 0, £ v = 

fit + fit) h 2 ~^ • ^ e em P nas i ze that there is only one co- 
herence length for the two order parameters when r — > 0. The 
London penetration depth is still given by Eq. (jTli with the 
order parameters given in Eqs. (46 1 and (47 I, Calculations of 
the order parameters near H C 2 and of the spatial correlation of 
(*¥ i(r)*¥ j(0)) when thermal fluctuations are present yield only 
one divergent coherence length consistently. Therefore, the 
superconductors are either type I or type II. 

Since it was discussed that vortex attraction originates from 
the overlapping of normal vortex cores,[ll | let us check the 
temperature dependence of the structure of a single vortex 
ranging from to T c , For simplicity, we take aff) - 
ff,(0)(l - T/Td) and assume yS, and y are T-independent. The 
parameters at T — are the same as those in the previous 
sections. The vortex structure is shown in Fig. [8] for several 
typical temperatures. It becomes clear that the sizes of vortex 
cores for Wi and v r , 2 get closer to each other as t approaches 0. 
Analytical calculations on the structure of the nonlinear vortex 
core also reveal an identical size for the vortex cores associ- 
ated with the two components close to r c .(36| Therefore, the 
interaction between vortices is purely repulsive in the present 
case. 

This result is understandable since the inter-band coupling 
is dominant at r <K 1, and thus the superconductivity in the 
two bands is strongly correlated and described by a unique co- 
herence length. Nevertheless, if the inter-band coupling y(T) 
and aj(T) all vanish at a single temperature T c as discussed in 
the previous sections, there will be two divergent lengths and 
thus different core sizes for the two components. This per- 
mits vortices to exhibit the peculiar interaction profile shown 
in Fig. B] Therefore, we find that the non-monotonic inter- 
vortex interaction, and thus the peculiar phenomena discussed 
above, only appear when the inter-band coupling and a^T) all 
vanish at T c within the framework of the GL theory. 



VII. CONCLUSION 

We have investigated the interaction between vortices in 
two-band superconductors. When the size of the magnetic 
flux of a vortex lies between the sizes of the two normal 
cores associated with the two condensates, the interaction be- 
tween vortices is attractive at large separation and repulsive at 
small separation. This was demonstrated clearly by the time- 
dependent Ginzburg-Landau calculations on two vortices in a 
superconducting disk. The equilibrium distance between two 
vortices is independent of the size of the disk for simulations, 
which clearly shows the existence of energy minimum in the 
interaction potential. The full dependence of interaction en- 
ergy on the vortex separation is obtained by the variational 
calculations. The two methods give the same estimate on the 
vortex separation with minimal free energy. 
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We have studied stable vortex configurations for a large 
number of vortices by Langevin dynamics adopting the novel 
distance dependence of vortex interaction. Circular vortex 
clusters coexisting with Meissner phase are observed for small 
and intermediate vortex densities. The transition from Meiss- 
ner state into vortex states is therefore of first order associated 
with a sharp increase of magnetization. The superconductiv- 
ity associated with the triangular vortex lattice is suppressed 
by a strong magnetic field in the same way of a type II super- 
conductor. Therefore, the magnetic behavior of these super- 
conductors as summarized in the mean-field phase diagrams 
interpolates those of the type I and type II superconductors. 
In most temperature regions except very close to the critical 
point, thermal fluctuations are weak for small Ginzburg num- 
ber, the same as single-band superconductors. 

The above interesting phenomena take place provided, first, 
the bare mean-field critical points of the two bands coincide 
with each other and the inter-band coupling vanishes at the 
same temperature, which permit two divergent core sizes as- 
sociated with the two condensates even close to the critical 
point; second, the parameters in the Ginzburg-Landau free en- 



ergy are appropriate to achieve the relations among the two 
coherence lengths and the penetration depth. 

When two condensates have different bare transition tem- 
peratures and the inter-band coupling is finite, the supercon- 
ductivity is induced by the inter-band scattering when temper- 
ature is sufficiently close to the critical point of the composite 
system. In this case, superconductivity in the two bands cou- 
ples strongly, and only one divergent coherence length exists, 
and thus the superconductors are either type I or type II. 
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